
clc;

disp('Dickson')
display('Voltage vectors')
% [in C5 C4 C3 C2 C1 out]
A1=[1 -1 0 0 0 0 -1; 0 0 1 -1 0 0 -1; 0 0 0 0 1 -1 -1; 1 0 0 0 0 0 0];
A2=[0 1 -1 0 0 0 -1; 0 0 0 1 -1 0 -1; 0 0 0 0 0 1 -1; 1 0 0 0 0 0 0];

A3=[1 0 0 0 0 0 0; 0 1 0 0 0 0 0;...
    0 0 1 -1 0 0 -1; 0 0 0 0 1 -1 -1];
A4=[1 0 0 0 0 0 0; 0 0 0 0 0 1 0;...
    0 1 -1 0 0 0 -1; 0 0 0 1 -1 0 -1];
NullA1_0=null(A1)
NullA1=NullA1_0(1:end-1,1:end);   % remove the last entry (corresponding to Vout)
NullA2_0=null(A2)
NullA2=NullA2_0(1:end-1,1:end);   % remove the last entry
NullA3_0=null(A3)
NullA3=NullA3_0(1:end-1,1:end);
NullA4_0=null(A4)
NullA4=NullA4_0(1:end-1,1:end);


% w vectors in the paper
dv1=NullA1_0;
dv2=NullA2_0;
dv3=NullA3_0;
dv4=NullA4_0;
 
dvc1=dv1(2:end-1,1:end)
dvc2=dv2(2:end-1,1:end)
dvc3=dv3(2:end-1,1:end)
dvc4=dv4(2:end-1,1:end)
%% Charge vector
display('Charge vectors')
% [in C5 C4 C3 C2 C1 out]
B1=[-1 -1 0 0 0 0 0; 0 0 -1 -1 0 0 0; 0 0 0 0 -1 -1 0; ...
    0 1 0 1 0 1 -1; 1 0 1 0 1 0 1];  % Phase 1a
B2=[-1 0 0 0 0 0 0; 0 -1 -1 0 0 0 0; 0 0 0 -1 -1 0 0; ...
    0 0 1 0 1 -1 -1; 1 1 0 1 0 1 1];  % Phase 2a
B3=[-1 -1 0 0 0 0 0; 0 0 -1 -1 0 0 0; 0 0 0 0 -1 -1 0; ...
    0 1 0 1 0 1 -1; 1 0 1 0 1 0 1; 1 0 0 0 0 0 0];  % Phase 1b
B4=[-1 0 0 0 0 0 0; 0 -1 -1 0 0 0 0; 0 0 0 -1 -1 0 0; ...
    0 0 1 0 1 -1 -1; 1 1 0 1 0 1 1; 0 0 0 0 0 1 0];  % Phase 2b

NullB1_0=null(B1)
NullB1=NullB1_0(1:end-1,1:end);     % again remove the last element
NullB2_0=null(B2)
NullB2=NullB2_0(1:end-1,1:end);
NullB3_0=null(B3)
NullB3=NullB3_0(1:end-1,1:end);
NullB4_0=null(B4)
NullB4=NullB4_0(1:end-1,1:end);

q1=NullB1_0;
q2=NullB2_0;
q3=NullB3_0;
q4=NullB4_0;

qc1=q1(2:end-1,1:end)
qc2=q2(2:end-1,1:end)
qc3=q3(2:end-1,1:end)
qc4=q4(2:end-1,1:end)


%% capacitor values
c=[1 1 1 1 1]';

%% Find common nullspace
dvc1=dvc1.*[c c c]
dvc2=dvc2.*[c c c]
dvc3=dvc3.*[c c c]
dvc4=dvc4.*[c c c]

NullQ1=null([dvc1 -qc1])
qc1=qc1*NullQ1(4:6)
NullQ2=null([dvc2 -qc2])
qc2=qc2*NullQ2(4:6)
NullQ3=null([dvc3 -qc3])
qc3=qc3*NullQ3(4:5)
NullQ4=null([dvc4 -qc4])
qc4=qc4*NullQ4(4:5)


%% Steady-state equation
NullQ=null([qc1 qc2 qc3 qc4])
% break
%% Final values

q1a=q1*NullQ1(4:6)*NullQ(1)
q2a=q2*NullQ2(4:6)*NullQ(2)
q1b=q3*NullQ3(4:5)*NullQ(3)
q2b=q4*NullQ4(4:5)*NullQ(4)

% normalized charge vectors
q=q1a+q2a+q1b+q2b;
qout=q(end);
q1a=q1a/qout
q2a=q2a/qout
q1b=q1b/qout
q2b=q2b/qout
